Thalamocortical Connectivity Development along the Sensorimotor-Association Axis: Sensitivity Analyses
Read in Data for Developmental Characterization
Glasser parcel names
#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv")
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) S-A axis
S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = c("orig_parcelname"), sort = F)Spatial permutation (spin) test parcel rotation matrix
perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regionsDataset-specific tract lists
#generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv")
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") Sensitivity variable dfs
#Surface area
area.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/cortical_anatomy/individual/PNC/PNC_surfacearea_finalsample.csv")
area.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/cortical_anatomy/individual/HCPD/HCPD_surfacearea_finalsample.csv")
#SNR
snr.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_SNR_finalsample.csv")
snr.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_SNR_finalsample.csv")
#Overlap sensitivity
sensitivity.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_overlapsensitivity_finalsample.csv")
sensitivity.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_overlapsensitivity_finalsample.csv")Sensitivity analysis results
setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/sensitivity_results/") #gam output dir, from gam_functions/fit_sensitivityGams.R
#list all .csv files in the environment results directory
files <- list.files(getwd(), pattern = '.csv', ignore.case = T, full.names = F)
for(i in 1:length(files)){
filename <- files[i]
Rname <- gsub('^.{12}|.{4}$', '', filename)
x <- read.csv(filename, header = TRUE)
assign(Rname, x)
}
rm(x)Association Between Tract-wise Sensitivity Variables and Main Analysis Maturational Timing
Main analysis thalamocortical development results
gameffects_FA_glasser_pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_gameffects_FA_glasser_pnc.csv")
gameffects_FA_glasser_hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_gameffects_FA_glasser_pnc.csv")Sensitivity variables in thalamocortical connections
tract_average_measures <- function(input.df, average.measure, final.tractlist){
input.df <- input.df %>% select(contains("autotrack")) #ensure we only have tract columns
tract.averages <- colMeans(input.df, na.rm = TRUE) %>% as.data.frame %>% set_names(sprintf("mean_%s", average.measure)) %>% mutate(tract = row.names(.))
tract.averages <- tract.averages[tract.averages$tract %in% final.tractlist$tract,]
tract.averages <- merge(tract.averages, S.A.axis, by = "tract", sort = F)
return(tract.averages)
}#surface area
regional.area.pnc <- tract_average_measures(input.df = area.glasser.pnc, average.measure = "area", final.tractlist = tractlist.pnc)
regional.area.hcpd <- tract_average_measures(input.df = area.glasser.hcpd, average.measure = "area", final.tractlist = tractlist.hcpd)
#connection SNR
tract.snr.pnc <- tract_average_measures(input.df = snr.glasser.pnc, average.measure = "SNR", final.tractlist = tractlist.pnc)
tract.snr.hcpd <- tract_average_measures(input.df = snr.glasser.hcpd, average.measure = "SNR", final.tractlist = tractlist.hcpd)
#atlas overlap sensitivity
tract.overlap.pnc <- tract_average_measures(input.df = sensitivity.glasser.pnc, average.measure = "sensitivity", final.tractlist = tractlist.pnc)
tract.overlap.hcpd <- tract_average_measures(input.df = sensitivity.glasser.hcpd, average.measure = "sensitivity", final.tractlist = tractlist.hcpd)#Combine development results and sensitivity measures
df.list.pnc <- list(gameffects_FA_glasser_pnc, regional.area.pnc, tract.snr.pnc, tract.overlap.pnc) #dfs to merge
gameffects_FA_glasser_pnc <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list.pnc)
gameffects_FA_glasser_pnc <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = join_by(orig_parcelname, label, tract))
df.list.hcpd <- list(gameffects_FA_glasser_hcpd, regional.area.hcpd, tract.snr.hcpd, tract.overlap.hcpd) #dfs to merge
gameffects_FA_glasser_hcpd <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list.hcpd)
gameffects_FA_glasser_hcpd <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = join_by(orig_parcelname, label, tract))Associations
surface area
cor.test(gameffects_FA_glasser_pnc$smooth.increase.offset, gameffects_FA_glasser_pnc$mean_area, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_pnc$smooth.increase.offset and gameffects_FA_glasser_pnc$mean_area
## S = 1273292, p-value = 0.5268
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.04500699
perm.sphere.p(x = gameffects_FA_glasser_pnc$smooth.increase.offset, y = gameffects_FA_glasser_pnc$mean_area, perm.id = perm.id.full, corr.type = "spearman")## [1] 0.35995
cor.test(gameffects_FA_glasser_hcpd$smooth.increase.offset, gameffects_FA_glasser_hcpd$mean_area, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hcpd$smooth.increase.offset and gameffects_FA_glasser_hcpd$mean_area
## S = 1238635, p-value = 0.5516
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.04256293
perm.sphere.p(x = gameffects_FA_glasser_hcpd$smooth.increase.offset, y = gameffects_FA_glasser_hcpd$mean_area, perm.id = perm.id.full, corr.type = "spearman")## [1] 0.37435
connection SNR
cor.test(gameffects_FA_glasser_pnc$smooth.increase.offset, gameffects_FA_glasser_pnc$mean_SNR, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_pnc$smooth.increase.offset and gameffects_FA_glasser_pnc$mean_SNR
## S = 1365701, p-value = 0.7327
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0243017
perm.sphere.p(x = gameffects_FA_glasser_pnc$smooth.increase.offset, y = gameffects_FA_glasser_pnc$mean_SNR, perm.id = perm.id.full, corr.type = "spearman")## [1] 0.44055
cor.test(gameffects_FA_glasser_hcpd$smooth.increase.offset, gameffects_FA_glasser_hcpd$mean_SNR, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hcpd$smooth.increase.offset and gameffects_FA_glasser_hcpd$mean_SNR
## S = 1063898, p-value = 0.01229
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1776307
perm.sphere.p(x = gameffects_FA_glasser_hcpd$smooth.increase.offset, y = gameffects_FA_glasser_hcpd$mean_SNR, perm.id = perm.id.full, corr.type = "spearman")## [1] 0.1184
atlas overlap sensitivity
cor.test(gameffects_FA_glasser_pnc$smooth.increase.offset, gameffects_FA_glasser_pnc$mean_sensitivity, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_pnc$smooth.increase.offset and gameffects_FA_glasser_pnc$mean_sensitivity
## S = 1135068, p-value = 0.03563
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1486775
perm.sphere.p(x = gameffects_FA_glasser_pnc$smooth.increase.offset, y = gameffects_FA_glasser_pnc$mean_sensitivity, perm.id = perm.id.full, corr.type = "spearman")## [1] 0.14175
cor.test(gameffects_FA_glasser_hcpd$smooth.increase.offset, gameffects_FA_glasser_hcpd$mean_sensitivity, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hcpd$smooth.increase.offset and gameffects_FA_glasser_hcpd$mean_sensitivity
## S = 1065814, p-value = 0.01305
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.1761501
perm.sphere.p(x = gameffects_FA_glasser_hcpd$smooth.increase.offset, y = gameffects_FA_glasser_hcpd$mean_sensitivity, perm.id = perm.id.full, corr.type = "spearman")## [1] 0.10545
Developmental Sensitivity Analyses
Number of significant developmental effects
Function to calculate FDR-corrected significant age effects
sig.effects <- function(dataset, sensitivity_measure){
ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset))
ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05
sigeffect.totaln <- sum(ageeffects.df$significant)
sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2)
print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes when controlling for %s", dataset, sigeffect.totaln, sigeffect.percent, sensitivity_measure))
}Sensitivity analysis: surface area
## [1] "In pnc, 202 thalamocortical connections (90.58 percent) show significant age-related changes when controlling for area"
## [1] "In hcpd, 180 thalamocortical connections (78.26 percent) show significant age-related changes when controlling for area"
Sensitivity analysis: connection SNR
## [1] "In pnc, 201 thalamocortical connections (90.13 percent) show significant age-related changes when controlling for SNR"
## [1] "In hcpd, 179 thalamocortical connections (77.83 percent) show significant age-related changes when controlling for SNR"
Sensitivity analysis: atlas overlap sensitivity (true positives)
## [1] "In pnc, 203 thalamocortical connections (91.03 percent) show significant age-related changes when controlling for sensitivity"
## [1] "In hcpd, 177 thalamocortical connections (76.96 percent) show significant age-related changes when controlling for sensitivity"
Cortex-wide thalamic connectivity developmental profiles
Function to plot tract-wise GAM smooths
developmental.smooths <- function(dataset, sensitivity_measure){
#get data for plotting
ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset))
smooths.df <- get(sprintf("smoothcentered_FA_%s_%s", sensitivity_measure, dataset))
smooths.df <- merge((smooths.df %>% filter(grepl("thalamus_R", tract))), (ageeffects.df %>% select(tract, GAM.smooth.partialR2)), by = "tract")
#match plotting params used in the main manuscript figures
if(dataset == "pnc"){
limit.low = 0
limit.high = 0.1
}
if(dataset == "hcpd"){
limit.low = -0.02
limit.high = 0.075
}
#plot tract-wise developmental trajectories
smooths.plot <- ggplot(smooths.df, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) +
geom_line(linewidth = .3, alpha = .8) +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(limit.low, limit.high), oob = squish) +
theme_minimal() +
labs(x="\nAge", y="Connection FA\n") +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))
ggsave(filename = sprintf("/cbica/projects/thalamocortical_development/figures/images/Figure_sensitivity_supplement/agesmooths_%s_%s.pdf", sensitivity_measure, dataset), device = "pdf", dpi = 500, width = 1.85, height = 1.65)
return(smooths.plot)
}Sensitivity analysis: surface area
Sensitivity analysis: connection SNR
Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis
Function to compute age of maturation correlation with the S-A axis and its spin-based significance
SA.alignment.statistics <- function(dataset, sensitivity_measure){
ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset))
ageeffects.df <- merge(ageeffects.df, (S.A.axis %>% select(tract, SA.axis)), by = "tract", sort = F)
spin.data <- left_join(spin.parcels, ageeffects.df, by = c("tract"))
corr.value <- cor.test(spin.data$SA.axis, spin.data$smooth.increase.offset, method = c("spearman"), exact = F)$estimate
pspin.value <- perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman")
results <- cbind(corr.value, pspin.value)
return(results)
}Function to plot thalamocortical connection age of maturation along the S-A axis
SA.alignment.plot <- function(dataset, sensitivity_measure){
ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset))
ageeffects.df <- merge(ageeffects.df, (S.A.axis %>% select(tract, SA.axis)), by = "tract", sort = F)
SA.plot <- ggplot(ageeffects.df, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
geom_point(size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
labs(x="\nS-A axis", y="Age of maturation (years)\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))
ggsave(filename = sprintf("/cbica/projects/thalamocortical_development/figures/images/Figure_sensitivity_supplement/SAaxis_alignment_%s_%s.pdf", sensitivity_measure, dataset), device = "pdf", dpi = 500, width = 1.85, height = 1.65)
return(SA.plot)
}Sensitivity analysis: surface area
## corr.value pspin.value
## rho 0.5114256 6e-04
## corr.value pspin.value
## rho 0.4913446 0.00185
Sensitivity analysis: connection SNR
## corr.value pspin.value
## rho 0.5046296 0.00085
## corr.value pspin.value
## rho 0.5129782 0.00215
Sensitivity analysis: atlas overlap sensitivity (true positives)
## corr.value pspin.value
## rho 0.5293327 0.00055
## corr.value pspin.value
## rho 0.5601407 0.0011